Exact results on the quench dynamics of the entanglement entropy in the toric code 
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We study quantum quenches in the two-dimensional Kitaev toric code model and compute exactly the time- 
dependent entanglement entropy of the non-equilibrium wave-function evolving from a paramagnetic initial 
state with the toric code Hamiltonian. We find that the area law survives at all times. Adding disorder to the 
toric code couplings makes the entanglement entropy per unit boundary length saturate to disorder-independent 
values at long times and in the thermodynamic limit. There are order-one corrections to the area law from 
the corners in the subsystem boundary but the topological entropy remains zero at all times. We argue that 
breaking the integrability with a small magnetic field could change the area law to a volume scaling as expected 
of thermalized states but is not sufficient for forming topological entanglement due to the presence of an excess 
energy and a finite density of defects. 



I. INTRODUCTION 

Cold atomic systems have provided an ideal playground for 
the experimental studies of the unitary evolution of thermally 
isolated quantum systems.'"-' Apart from the experimental in- 
terest, studies of quantum evolution have shed light on some 
fundamental questions such as thermalization4i^ 

A widely-studied scenario is the unitary evolution of a 
quantum state when the Hamiltonian is suddenly changed - 
a "quantum quench"^— which has been the subject of many 
recent studies (see Refs. Ill2l - [l5ll for a few examples). A chal- 
lenging question, which has been addressed theoretically in 
one dimensional systems, is the evolution of the entanglement 
entropy following the quench.'^''' More generally, such stud- 
ies are concerned with the fate of an out-of-equilibrium quan- 
tum state with a Hamiltonian that changes with some time- 
dependent protocol. 

The case where the coupling constants are changed across 
a quantum phase transition is of special interest. In recent 
experiments with spinor condensates, for example, the forma- 
tion of ferromagnetic domains was observed in a quench to a 
Hamiltonian with a symmetry-breaking ferromagnetic ground 
state. ^ Even though the interactions after the quench favor a 
certain type of order, it is not at all obvious that the order 
should emerge out of equilibrium. Since symmetry-breaking 
order has a local order parameter, the non-equilibrium state 
can support the formation of ordered domains. But what about 
quenches across a topological phase transition? Would any of 
the topological characteristics of the equilibrium phase show 
up in the non-equilibrium wave-function? 

Another interesting aspect of the non-equilibrium dynam- 
ics is the issue of thermalization and whether a large system 
can act as an effective heat bath for its subsystems. If thermal- 
ization occurs, due to the presence of some excess energy, the 
non-equilibrium state at long times should be in some sense 
close to a thermal state at a finite temperature, which implies 
volume scaling for the entanglement entropy. 

The qualitative discussions above motivate a detailed anal- 
ysis of a quantum quench across a topological phase transition 
into the topological phase. We study the specific case of the 
Kitaev toric code model'** with spin-polarized initial states. A 
related quantum quench problem with topologically ordered 
initial states was recently studied numerically in Ref. II19I1 . We 



focus on the time evolution of the entanglement entropy be- 
cause it is intimately related to topological order through the 
topological entanglement entropy i22i21 Entanglement entropy 
is also related to thermalization since, as mentioned above, 
volume scaling of the entropy as opposed to the area law is 
generically expected of thermal states. 

We compute exactly the bipartite entanglement entropy as 
a function of time. This is an example where exact calculation 
can be done in 2h-1D (thus far, analytical results for the time- 
evolution of entanglement had focused on 1h- ID systems). We 
show that the entanglement entropy respects the area law at 
all times and if one adds disorder to the coupling constants, 
the recurrences or revivals are suppressed at long times and 
the entanglement entropy per unit "area" (boundary length) P 
equals 

S/P = ln(4/e) 

in the thermodynamic limit when the initial polarization is in 
the x or z direction. 

We find the dependence of this saturation value on the ini- 
tial state for an arbitrary du^ection of the initial magnetiza- 
tion. A new ingredient that comes in when considering these 
initial states is order-one corrections to the entanglement en- 
tropy from the convex and concave corners in the subsystem 
boundary. We find that the general form of the generated en- 
tanglement entropy is given by 

5(t) = /(i)P + />(i)C>+/<(i)C< 

where P is the perimeter of the subsystem and C> (C<) the 
number of convex (concave) corners in its boundary. The 
functions of time /, /> and /< are independent of the sub- 
system geometry. From the form of the entanglement entropy 
above, we find that the topological entropy vanishes at all 
times. 

We study the effects of breaking the integrability in a per- 
turbative approach. We identify the mechanism by which an 
integrability-breaking magnetic field can change the area law 
to the volume law as expected of thermalized states. The phys- 
ical picture is that the area law survives at low orders in an 
expansion in the integrability-breaking field. Going to higher 
orders however increases the thickness of a ring-like region 
around the subsystem boundary that contributes to the entan- 
glement entropy. Eventually at high enough orders, the ring 



will collapse onto itself making volume scaling possible. In 
the perturbative treatment, we explicitly find that the topo- 
logical entanglement entropy remains zero at low orders. We 
expect that the excess energy and the finite density of defects 
should prevent the formation of any topological order in a sud- 
den quench even if one breaks the integrability. 

The paper is organized as follows. In section|II]we discuss 
quantum quenches for spin-polarized initial states in the x or 
z direction. We calculate the entanglement entropy analyti- 
cally and discuss the effects of disorder in the coupling con- 
stants. In section [nil we study the role of the sharp corners in 
the subsystem boundary in the generation of entanglement by 
exactly calculating the second order Renyi entropy for an ar- 
bitrary direction of the initial magnetization. We then discuss 
the effects of breaking the integrability, which we treat pertur- 
batively, in section |IV] We conclude the paper in section [Vl 



II. INTEGRABLE QUENCH: INITIAL MAGNETIZATION 
IN THE X OR z DIRECTION 

Consider the toric code Hamiltonian in a field22.— 



if = - ^ XpBp - ^ ^isAs - hy^/i 



(1) 



The Hilbert space of the above Hamiltonian is given by spins 
living on the bonds of a square lattice on a torus and the sums 
are over all plaquettes, stars and bonds respectively. The pla- 
quette and star operators are defined as follows 



Br, 






A.. 






where as are the Pauli matrices. For the usual Kitaev model 
the coupling constants A and /i are independent of the plaque- 
ttes and the stars but here we keep the p and s indeces to allow 
for disorder in the couplings. Note that for the model in zero 
magnetic field (h = 0), the ground state is independent of the 
actual values of the coupling constants A and /z as long as they 
are all positive. 

One can study the dynamics of this model by changing the 
coefficients A, /i and h in time with some protocol. In this 
section, we consider a quantum quench where the system is 
initially in a fully magnetized state and then evolves with the 
integrable toric code Hamiltonian, namely, we are initially in 
the ground state for h ^ and Ap = /Zs = for all stars 
and plaquettes. We then turn off the field h at t — (which 
does not change the quantum state) and turn on the toric code 
part of the Hamiltonian as follows: A == X{t) and /i = /x(t) 
(A(0) = fi{0) = 0). 

Consider for now, the case where n — x {n ■ a — a^). 
The ground state of the usual model (without disorder in the 
couplings) exhibits a topological phase transition: for h > he 
the system is a paramagnet and for h < he, it is topologi- 
cally ordered. This topological phase is robust against weak 
disorder in the couplings. Let us call the initial state |17). In 
the case where the initial magnetization is in the x direction 



irix) = |11...) in the |crf erf ...) basis. The initial state is a ten- 
sor product of single spin states and thus has no entanglement. 
In the paramagnetic initial state we have (a^) = 1 for all 
spin operators. Before addressing the question of the entan- 
glement entropy evolution, let us consider the simple problem 
of what happens to the <j^ expectation value. As a function 
of time, we have (o'^(t)) = {ilx\W(J^U\ilx) where U is the 
evolution operator given by 



U = exp 



Z^ 



\p{t')dt'Bp + i 



•E 



^is{t')dt'A, 



(2) 
Since all star and plaquettes commute with one another, we 
have dropped the time ordering operator. If we consider a sud- 
den quench, we have Jg Xp{t')dt' = Xpt and J^ iis{t')dt' = 
fist- Let us focus on sudden quenches for a more compact 
notation. At the end of the calculation we can replace Xt and 
/it with the corresponding integrals over time-dependent cou- 
plings for a more general protocol. 

The expectation value of a^ , , a spin between two nearest- 
neighbor plaquettes p and p' is then given by 



{a;p,{t))^{nx\e 



itXj, 



-it\„iB. 



p' "p' cr^p, e^*^"-^" e** V -Bp' |f7^) . 



Since Bz, — 1, we have 



i iL/^j} -D T) 



cos(Apt) ± i sin{Xpt)Bp 



We have B„ „> a. 



p-.p "pp' 



-Upp,-Dp,p'. 



Therefore we get 



{a^^, (t)) = cos (2Apt) cos (2Ap't) . 

In the completely clean case (a theoretical abstraction), the 
average of (0-^ , {t)) over time is equal to i but any amount of 
disorder will eventually lead to a state which unlike the initial 
paramagnet has a vanishing time-averaged expectation value 
for a'^ . 

We next consider the main problem of this section, namely 
the evolution of the entanglement entropy. The density matrix 
of the system can be written as p{t) = U\il){il\U'' where U 
is the time evolution operator and \il) is a generic initial state. 
By inserting the resolution of identity twice, we can write. 



pit)^Y.('^\U\n){n\U^a)\a){a\ 



(3) 



We would like to calculate the entanglement entropy be- 
tween two subsystems A and B. Any star or plaquette opera- 
tor acts on four spins which may all belong to A, all belong to 
B or some belong to A and some to B. In this last case, the 
star or plaquette operator must lie at the boundary of the two 
subsystems. We can then write 



U = UaUbUo 



(4) 



where Ua and Ub only act on spins in A and B respec- 
tively and Ug only depends on star and plaquette operators 
at the boundary of the two subsystems. The three opera- 
tors Ua, Ub and Ug commute with one another. Each ba- 
sis state can be written as the tensor product of states for 



spins in A and B, \a) = laAPs)- We now integrate out the 
spins in B to find the reduced density matrix for subsystem A, 
PA{t) = E/3b(/3bIp(*)I/3b) and obtain. 






{aA/3B\U\n){n\U^aApB)\aA){ 



<7a\ 



(5) 

Notice that since the trace of an operator is independent of the 
basis, we can replace |/3b) by Ub\(3b) in the above expres- 
sion. This amounts to calculating the trace in a rotated basis. 
We then substitute Eq. ^ into Eq. ^ above and obtain 

PAit)= Y. iaApB\UAUa\^){mlul\aApB)\aA){aA\ 

The part of the Hamiltonian that only acts on B spins has no 
effect on the time evolution of the reduced density matrix for 
subsystem A. In order to calculate the entanglement entropy 
between the two subsystems we need to find tr(p^). We give 
a replica index to a, /3 and a and drop the subsystem index for 
brevity. States \ai) are in subsystem A and \/3i) in B. Using 
{ai\ai+i) = Sai,ai+i we obtain. 



tr(pl)= Y ll'^»^^^\UAUa\n){n\uiul\a,+,p,) 

{a,},{ft}i=i 

(6) 
with an+i = «!. Once again because traces are basis- 
independent we can replace | a; ) by Ua | a^ ) in the above equa- 
tion and we obtain an expression for tr(/9^) that only depends 
on Uq and not the whole U. Let us define the following com- 
muting operators 

U± = e±^^peo, KB,t^ ^± ^ gizE.ea, /^=^=* (7) 

where dp (dg) is the set of all boundary plaquettes (stars) 
which act on both A and B spins. We then obtain 



flips on \Vt,x)- Let us consider then Ising spin variables living 
at the center of dp plaquettes. A plaquette is flipped for 9 = 
— 1 and not flipped for 6 = 1. Notice that, on a torus, flipping 
a certain set of plaquettes {9} and flipping all other plaquettes 
that do not belong to {9} have the same effect on the spins, 
i.e. they describe the same element in the group of plaquette 
flips. However, since we have eliminated Ua and Ub from 
the formulation of the problem by a basis rotation, we can 
uniquely specify a state with nonvanishing {al3\Up\flx) by 
a set of flipped boundary plaquettes. Let the state ]ai(3i) be 
labeled by {9p} forp £ dp. If the partitions are not too small, 
it is not possible for a set of boundary plaquette flips acting on 
\flx) to create a state with the same B spins and different A 
spins or vice versa. Therefore \a2P1) = |cii/3i) = \{9p}) and 
we obtain using Eq. ([8|l 



tripl) 



E 

{8p},p€dp 



{{9p}\u+\nx) {nx\Up\{9p}) 



(10) 

Using Eq. (|9]l we find that {{9p}\Up\i^x) is a product of terms 
that are equal to cos Xpt for each unflipped plaquette with 
9p = 1 and i sin Apt for each flipped plaquette with 9p = —1. 
We can then write 

{{9p}\U+\nx) ^ n (^^ cos Apt + zi^ sin Apt 

(11) 
Switching the order of the sum and product after plugging 
Eq. ( fTTT i into Eq. ( fTOl i. we obtain 



p£dp 9p 



We have used the fact that for 9 = ±1, ^ = {^) . 

Summing over 9p = ±.1 for each boundary plaquette sim- 
ply gives 



{a,},{/3.}»=l 

(8) 
with a„+i = ai. We now focus on the simple case where 
initially all spins are aligned in the x direction. In this case, 
\Vt) = \Vtx) is an eigenstate of Us and we can therefore sub- 
stitute U+\Q.x){^x\Us = \^x){^x\ in Eq. (O above. Also 
I a) and |/3) are chosen in the a^ basis so the action of Bp on 
a basis state is to flip the four spins around the plaquette p. 
For concreteness, we assume that the partition B contains all 
spins which lie inside a simply connected M x N rectangular 
region. We then have 2(Af + N) — 4 boundary plaquettes. 
Let us consider the action of U^ = Yip^f^ e^p^^* on a state 



|q;/3). Since Bp = 1, we have 



i-pedp 



U^ = Yl (cos Apt ± iBp sin Apt) . (9) 

p£dp 

Notice that {al3\U+\Clx) ^ only if the state the state \a/3) 
is obtained by acting with a certain number of dp plaquette 



i^ipA) = n [(cosApt)2" + (sinApt)^"] . (12) 

p€dp 

By analytic continuation of n, we write the entanglement en- 
tropy as 

d 

S = -tr {pA log{pA)) = lim 7^tr(p^) 
n^i dn 

and we finally obtain, 

S = - Yj [^°^^ Apt In (cos^ Apt) + sin^ Apt In (sin^ Apt)] . 
pedp 

(13) 
For the clean Kitaev model, the entanglement entropy is an 
oscillatory function of time whose maximum is proportional 
to the number of boundary plaquettes or the perimeter of the 
subsystem. For this initial condition, the star term in the toric 
code Hamiltonian has no effect on the evolution of the entan- 
glement entropy because the initial state is an eigenstate of 
that term. Similarly if initially all spins are in the z direction. 



the plaquette term will have no effect. The completely clean 
system will never reach a steady state. 

If there is some disorder in the coupling constants however, 
the entanglement entropy per unit boundary length will satu- 
rate to a steady state value in the thermodynamic limit. Sup- 
pose the number of boundary plaquettes P is very large. As 
t — > oo, due to the presence of some disorder in couplings Ap, 
the phases Xpt are equally likely to take any value between 
and 277 modulo 27r. Therefore, taking the thermodynamic 
limit first we can write S{{4>i}) — — X]i=i C{4'i) with 

(((/)) = cos^ 01n (cos^ (f)) + sin^ cfiln (sin^ (f) 

where (j)i is a random variable between and 2tt drawn from 
a uniform distribution. The average entanglement entropy is 
then given by 

r27T 

d(j)C{(f))=P\n{4:/e) (14) 



D 



27r 



We have (5 
of disorder 



5)-^ ex P which implies that for any realization 



I p p 



\-0{^=) 



P 



indicating that in the thermodynamic limit, the entanglement 
per unit area length saturates to a steady-state value ln(4/e). 
Notice that these results, namely the survival of area law at 
all times and the oscillatory behavior in the entanglement en- 
tropy in the absence of disorder are generic features of inte- 
grable quenches where the integrals of motion are operators 
with local support. The entanglement entropy per unit bound- 
ary length relaxes in the thermodynamic limit to a steady state 
value in the presence of disorder in the coupling constants. 

Because the system can reach a steady state, it is natural to 
ask if the steady state of a subsystem has the properties of a 
thermalized state described by a Gibbs or a generalized Gibbs 
ensemble. ^^ For such thermalized states however, we expect 
the entanglement entropy of a small subsystem to scale as its 
volume in the limit of long times. The survival of the area law 
at all times which follows trivially from having local integrals 
of motion indicates the absence of thermalization even to a 
generalized Gibbs ensemble. 
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FIG. 1. Subsystem B consists of spins inside the painted region, 
while subsystem A includes spins outside and on the boundary of 
this region. The convex (concave) comers for subsystem B are rep- 
resented by squares (circles). 



pointsi^— Here we find that corners can also make sublead- 
ing contributions to the dynamical entanglement. We find that 
when subsystem B consists of spins inside a region with C> 
convex and C< concave corners as seen in Fig. [T] the entan- 
glement entropy after a quantum quench behaves as 



5-/P + />C> + /<C< 



(15) 



where P is the perimeter of the subsystem and /, /> and /< 
are functions of time that do not depend on the subsystem 
geometry. For a subsystem with a mostly smooth boundary 
and only a few sharp corners, the leading term is the fP but 
for a subsystem with a corrugated boundary, the entanglement 
entropy can be different in the thermodynamic limit. 

Notice that the number of concave and convex corners are 
not independent. Suppose subsystem B consists of spins in- 
side of an area with m connected regions, each with a cer- 
tain number of handles such that the total number handles is 
g. As an example in Fig. [T] we have C> = C< = 11 and 
TO = g = 2. We argue that 



C, 



C<+4(to-(?). 



This is because for a rectangular simply connected region, we 
have C> = 4 and C< — and transforming the shape by 
adding or subtracting plaquettes on the square lattice does not 
change C> — C< . Therefore each simply connected region 
makes a contribution of +4 to C> — C< and by a similar 
argument, each hole contributes —4. 



III. SUBSYSTEM CORNERS AND THE DYNAMICAL 
ENTANGLEMENT ENTROPY 



A. Initial magnetization in the y direction 



Let us now consider a more generic case where the initial 
density matrix is not diagonal in either cr^ or cr^ basis and both 
terms in the Kitaev Hamiltonian affect the entanglement dy- 
namics. We find that even though in the thermodynamic limit 
the area law still applies, the physics of the out of equilibrium 
entanglement is much richer. We find that here there are order- 
one corrections to the entanglement entropy from the concave 
and convex corners of the subsystem. 

Subleading corrections to the entanglement entropy aris- 
ing form the sharp corners in the boundary were found be- 
fore for static entanglement entropy in the logarithmic correc- 
tions to the area law in a class of conformal quantum critical 



To see how the contribution from sharp corners comes in, 
consider an initial state with all the spins in the positive y di- 
rection. Let us work in the the (7^ basis with cr*'|±) — ±|±), 
(T^|±) — |=f) and cr^|±) = T*|t)- The action of a star or 
plaquette operator in this basis is then to flip the correspond- 
ing four (jy spins modulo a phase factor of ±1 from Bp. We 
will see later that this sign has no effect on the entanglement 
entropy. 

We assume for simplicity that there is a minimum distance 
of 3 times the lattice spacing between all corners. Consider a 
convex corner as in the right hand side of Fig.|2l The /? spins 
live in the dashed blue (color online) bonds and the a spins 



flupn 
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II 



III 



IV 



FIG. 2. A convex (concave) comer of the B subsystem is shown on 
the right (left) hand side. The dashed blue bonds belong to subsystem 
B and the solid black bonds to subsystem A. A comer boundary 
plaquette (star) is represented by a red square (cross). 



on the solid black bonds. Flipping the corner plaquette repre- 
sented by a red square has the same effect on /3 spins as flip- 
ping the two boundary stars represented by red crosses. These 
operations of course have a different effect on a spins. Be- 
cause we assumed a minimum distance between the corners, 
these corner star and plaquette operators uniquely correspond 
to one corner. Similarly as in the left hand side of Fig.|2] flip- 
ping the corner star at a concave corner has the same effect on 
the a spins as flipping the two corresponding plaquettes. We 
see from this discussion that the effects of stars and plaquette 
flips mix at the corners, which leads to a different entangle- 
ment generation than along a smooth boundary. 

The details of the derivation will be given in appendix [A] 
Let us write for now the result for the entanglement entropy 
for a subsystem B consisting of the bonds inside an A^ x M 
rectangular region. A boundary star or plaquette that is not 
located at a corner is refereed to as a regular star or plaquette. 
As seen in Fig. |2] at each convex corner, one plaquette and 
two stars are special. For our N x M subsystem, we have four 
convex corners, 2(A/+A^— 4) regular boundary plaquettes and 
2{M + N — 6) regular boundary stars and the entanglement 
entropy for a clean system after a sudden quench is given by 

S = 2{M + N ^ 4)sp + 2{M + N - 6)s, + As^ (16) 

where the contributions from regular plaquettes (Sp), regular 
stars (ss) and convex corners (sc) are given by 

Sp = (sin^ Xt In sin^ Xt + cos^ Xt In cos^ Xt) 
Ss = (sin^ fit In sin^ fit + cos^ fj,t In cos^ ^t) 
Sc = [ri In n + 2r2 In r2 + ra In ra] 



with 



r2 
r3 



sin^ Xt sin*^ ^t + cos^ Xt CDs'* /zi 



sin fit cos ^t 



sin^ Xt CDs'* lit + cos^ Xt sin^ fit. 



The above solution conforms to the generic form Eq. (fTsT i for 
P = 2{M -I- TV - 2), C> = 4, C< = 0, / = Sp + s^ and 
f> ^ Sc- 8ss - 4sp. 

In appendix lAl we give a derivation of the above result in 
the more general setting of an arbitrary subsystem geometry. 
If instead of a sudden quench we turn on the coupling con- 
stants with some protocol A(i) and /i(t), we just need to re- 
place At and fit by their integrals over time. 



FIG. 3. The partitions used for calculating the topological entangle- 
ment entropy. 

Using the partitions of Ref. '2V, shown in Fig. [3] the topo- 
logical entanglement entropy is given by 



'S'topo — Si — Su — S] 



III 



s, 



IV- 



(17) 



From the form of Eq. ( fTsT i, we can then see that although there 
are order-one corrections to the area law, the topological en- 
tanglement entropy remains exactly zero at all times indepen- 
dently of how the toric code Hamiltonian is turned on. 

It seems plausible that for any initial state, as long as the 
integrability is not broken (with a field h for example), evo- 
lution with the Kitaev Hamiltonian does not change the topo- 
logical entanglement entropy for any time dependence of the 
coupling constants. This is because this integrable evolution 
does not allow propagation of information. 



B. Arbitrary initial magnetization 

The von Neumann entropy is a special case (limit of n 
of the order n Renyi entropy defined as 



1) 



Snip a) = 



1 



1 -n 



tripl). 



Hereafter we refer to the second order (n = 2) Renyi entropy 
as the Renyi entropy. The Renyi entropy 5*2 (pa) — ^tr(p^) 
is simpler to calculate than the von Neumann entropy and is 
an equally valid measure of entanglement. It has been shown 
that the topological entanglement entropy calculated by any 
order Renyi entropy contains the same information in the abe- 
lean topological phases^'' such as the toric code. For an in- 
tegrable quench from a paramagnetic state with magnetiza- 
tion in an arbitrary direction, we exactly calculate the Renyi 
entropy 5*2 (/Oa)- We consider the clean case with position- 
independent coupling constants. With disorder in the cou- 
plings, the entropy saturates to the time-average of the clean 
case similarly to the simpler initial conditions considered be- 
fore. 

In the ax basis, an arbitrary paramagnetic state can be writ- 
ten as 



\n) 



-)) 



^N 



The calculation is rather involved and is done by a transfer 
matrix method explained in appendix IB] We find the same 
general features for an arbitrary initial magnetization as for 
the initial magnetization in the y direction. Let us summarize 
these features below. 




FIG. 4. The dependence of the saturation value of the generated 
Renyi entropy per unit boundary length on the two angles parame- 
terizing the magnetization direction of the initial state. 



The Renyi entropy can be written as a sum of contributions 
from the perimeter of the subsystem as well as the convex and 
concave corners. These contributions are oscillatory functions 
of time in the absence of disorder in the couplings which at 
long times, saturate to the values given by the time-average 
of the clean case if we have some disorder in the couplings. 
The saturation values are independent of the final coupling 
constants and the protocol used for turning them on. They do 
however depend on the initial state. 

For a large subsystem with only a few sharp corners, the 
leading term for the Renyi entropy is the first term in Eq. ( fTSl ). 
namely S « fP. The saturation value of the Renyi entropy 
per unit boundary length is plotted in Fig.|4]as a function of 
the 0/2 and $ angles which determine the initial state through 
c^ = cos 8/2 e'* and c^ = sin 6/2. Notice that this satura- 
tion value is a unique function of (Bp) = (sin Q cos $) and 



(^) = (cose)^ 



IV. BREAKING THE INTEGRABILITY: 

ENTANGLEMENT ENTROPY IN THE INTERACTION 

PICTURE 

To study the effects of breaking the integrability, we fo- 
cus on one simple case. Namely, we start with a state that 
has all spins aligned in the x direction. As opposed to the 
integrable case where the magnetic field is turned off before 
turning on the integrable toric code Hamiltonian, the the toric 
code Hamiltonian is turned on in the presence of a small mag- 
netic field h. We cannot perform exact analytical calculations 
in this case. We can however discuss the effects of the mag- 
netic field by calculating a formal expansion for the (Renyi) 
entanglement and topological entropies in powers of h. The 
expansion can only be considered as a perturbative treatment 
if the toric code coupling constants are quickly pushed above 
h and the time scales considered are not too large. 

Two main insights are derived from this analysis. First we 



find indications as to how breaking the integrability can even- 
tually lead to volume scaling in the entanglement entropy. The 
physical picture is that in the expansion in powers of h, the 
contribution to the entanglement entropy comes from ring- 
shaped areas around the subsystem boundary whose thick- 
ness grows at higher and higher orders in h. Eventually, at 
a high enough order, these correlated regions meet and the 
terms in the expansion no longer respect the area law. Even 
for a small h, these high order terms become important at 
long times making volume scaling of the entanglement en- 
tropy possible as expected of thermalized states. Secondly, 
we find that in this expansion in powers of h, the topologi- 
cal entanglement entropy remains zero at low orders until we 
reach an order proportional to the system size, suggesting that 
a small integrability-breaking term is not enough for generat- 
ing topological entanglement. 

A useful notion that we introduce in this section is the en- 
tanglement entropy in the interaction picture. We are used to 
the Shcrodinger, Heisenberg and interaction pictures for cal- 
culating the expectation values of observables, which are of 
course independent of the picture. If one defines the entangle- 
ment entropy as a property of the wave function alone, then 
since the wave function does not change in the Heisenberg 
picture, calculating the entanglement entropy in that picture 
leads to the nonsensical result that time evolution does not 
change the entanglement of a quantum state. An interesting 
alternative is to try and define the entanglement entropy en- 
tirely in terms of the correlation functions of appropriate op- 
erators. We do not pursue this path here. We argue instead 
that as explained below, the entanglement entropy calculated 
only from the wave function in the interaction picture can be 
a useful concept. 

Let us first consider a non-interacting Hamiltonian Ho such 
as a paramagnet. Consider a partitioning of the system into 
subsystems A and B. The evolution operator for Ho can be 
written as 

Uoit)=UoAit)UoB{t) 

where Ua (Ub) only acts on the degrees of freedom in A (B). 
Notice that acting with Uo (t) on any wave function does not 
change its entanglement properties. Now consider a Hamilto- 
nian 

H = Ho + V 

with the corresponding evolution operator ZY(t). 

The wave functions at time t in the Schrodinger and inter- 
action pictures are respectively given by 

\Mt)) = uit) iV'(o)), \Mt)) = u^it) u{t) |v;(o)). 

Now since acting with Uo or Ul on any wave function does not 
change its entanglement, the interaction picture wave function 
has the same entanglement entropy as the Schrodinger picture. 
Note however the calculation may be much simpler in the in- 
teraction picture. 

In the discussion above, we made use of the interaction pic- 
ture entanglement entropy because we had a Hamiltonian Ho 



whose evolution operator did not change the entanglement en- 
tropy. At the next level of complexity, where using the in- 
teraction picture wave function may be useful for calculating 
the entanglement entropy, the evolution operator of the unper- 
turbed Hamiltonian {Uq) may have another known property. 
For example it may only create entanglement entropy scaling 
with the subsystem area or it could be incapable of changing 
the topological entanglement. This is the case for the inte- 
grable toric code Hamiltonian. Now if we are interested in 
finding whether a perturbation can generate entanglement en- 
tropy with volume scaling or whether it can generate topolog- 
ical entanglement entropy, we can perform the calculation for 
the wave function \tp\{t)) in the interaction picture. 

Let us go back to our quench problem and write the wave 
function in the interaction picture with 



and 



V^hY^^^r 



Starting with the initial state |0) = \Q,x) with all the spins 
in the positive x direction, the wave-function at time t in the 
interaction picture is given by 

|V.i(t)) =e''/o'i*'^o(t')|^g(i)) ^ ijy{t)\Q) (18) 

where Uy {t) is given by 

Uv{t) = l-ih I dtiVih) 



for 



+ i-ihf / dti / dt2V{h)V{t2) + 

Jo Jo 



V{t)^he-^fodt'Ht')j:,B, K-^. 



^ 1 (,ifodt'Ht')i:^B^ 



The sums are over all bonds, plaquettes and dominos and 
the graphical notation represents a product of Pauli matrices. 
For example if the bonds in a plaquette are numbered clock- 
wise from 1 to 4 with the top bond being number 1 we have 



^1^2 3 4 



^1^2 3 4 
(^z'^y<^z^zi 



The closed form in Eq. ( fT9l ) is simply derived from the re- 
peated substitution of the following commutation relations 







r 


r^p, c 


= 2iP, 


p J 


1 


r^Bp, C + D 


= 4iP 


p 


J 





E^p'^ 



-4i{C + D) 



in the Baker-Hausdorff expansion of V 



V{t) 



i-ixty 

2! 



h i-iXt) 



T.Bp.T.'' 
E^p'E< 



and resumming the series. 

The Renyi entropy for the interaction picture wave function 
\ip\{t)) is given by 

S2{t) - -log(ai/3i|?7y|0) {Q\Ul\a2l3i) x 

{a2P2\Uv\Q) ml\aili2) (21) 

where summation over repeated indices is implied and as be- 
fore, the a and /? spins belong to the A and B subsystems 
respectively. First let us consider the first order term in the 
tr(p^) (the expression for the Renyi entropy without the log). 
We use the shorthand notation 



VV---V 



dti 



t„-i 



dtS{tl)---V{tn). 



For the simplicity of notation, we focus on the sudden quench 
and write J„ dt'X{t') as Xt, bearing in mind however that for 
a more general protocol, we can substitute the integral at the 
end of the calculation. 

We can explicitly calculate V{t) and the result is 



cos4At+l ^ sin4Ai ^ cos4At-l „ 
V{t) = C H z — P -\ z D 



2 2 

= c{t)C + p{t)P + d{t)D 

with the operators C, P and D defined as follows 

^ = E< 



(19) 



Ey z z z 

(z z-|-z y -\- ^ z + y z) 
^ |_ z z y z ' 

V 

^-ER 



(20) 



and similarly for the Hermitian conjugate. The first order term 
vanishes as seen below 



tr(p^) = 1 - 2i/i(0|V"|0) + 2i/i(0|V"t|o) 



l + 0(/i2 



The second order contribution to tr(p^) is given by 

tr(p^)L(,.) = h" [2 |(0/3|y|0)p + 2 |(a0|y|0)|2 - WW 

- (0|1/^|0)2 - 2 (0|yV^|0) - 2 (0|FVt|0) ] . 

Let us assume that the whole system is on a torus and the 
total number of bonds (spins) is Nh- The number of plaquettes 
Np = -1^ and the number of dominos in the system is Nd — 
Ni,. Let us introduce another shorthand notation for integrals 
involving the functions c, p and d, 



(/l,/2, ■ 



, /n) = / dti 

Jo 



dtnfl(ti)- ■■ fnitn) 
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so we have for example (c) = Jp dti c{ti) and {c,p) — 

Jp dti Jp^ dt2 c{ti) p{t2)- We denote the number of plaque- 
ttes (dominos) with all bonds in subsystem A by n^^ (ridA) 
and similarly by UpB ijidB) for subsystem B. We then have 

|(0^|y|0)|2 = 7V,2(^)2 _^ ^^^ X 16 (p)2 + n^Bid? 
|(aO|y|0)|2 = ^¥^2(^)2 ^ ^^^ X 16 (p)2 + n^Aidf 

(0|y|0) = iV,(c) 
(0|yy |0) - iVb2(c, c) + TVp X 16 (p,p) + NM d). 

The factors of 16 = — 4i x 4i come from the fact that each 
term in P is a sum of four operators that flip the spins around 
a plaquette. Noting that because c{ti) 0(^2) is obviously sym- 
metric under the exchange ti t2, we have (c, c) = (c)2/2 
and similarly for (p,p) and (d, d). We then obtain 



tr(p- 



= 1 - ft. (A'p - UpA - npB){p) 
- h^iNd - UdA - ndB){df, 
which using log(l + x) ^ x + 0{x^) gives 

S2 = h^iNp - Up A ~ ripB){pf 

+ h\Nd ~ UdA - ndB){df + 0(/i^). 



(22) 



This expression demonstrates the physical picture for building 
volume entanglement which we alluded to before. We observe 
from the above expression that at this order the contribution to 
the entropy is from a ring-shaped area of thickness two lattice 
spacings around the subsystem boundary. At higher orders, 
the thickness of this correlated region, which contributes to 
the entanglement entropy, grows until the ring closes in on 
itself, making volume scaling possible. 

Notice that in contrast to exact area law which is possible 
in a finite system, since the Renyi entropy is symmetric under 
exchanging the subsystems, volume scaling can only hold ex- 
actly for a finite subsystem in an infinite system and as a good 
approximation in the limit when a subsystem is much smaller 
than the whole system. 

It is a simple exercise in geometry to count np and Ud of 
the two subsystems for partitions shown in Fig. [3] and verify 
directly that no topological entanglement is generated at order 
h?. The details of this calculation will be given in appendixICl 
We will also show in that Appendix that at the next order with 
a non-vanishing contribution to the Renyi entropy (©(/i^)), 
the topological term will still vanish due to an exact cancella- 
tion. 

Similarly to the volume scaling, topological entanglement 
is non-perturbative and does not show up in an expansion in 
powers of h in the thermodynamic limit. For a finite system, 
we expect non-vanishing contributions to the topological en- 
tropy to appear at orders that scale with the system size. Let us 
comment here that on physical grounds, in the case of a sud- 
den quench, we do not expect any topological entanglement 
to form even at infinitely long times. A trend toward forming 
volume entanglement observed at low orders suggests that in 
the limit of large times, the generic thermalization behavior is 
expected. Although we cannot give a mathematical proof, the 
fragility of topological order to thermal fluctuation&i^Ml and 



the extensive excess energy in the system suggests that the 
non-equilibrium state, after a sudden quantum quench, will 
likely bear no signature of the equilibrium topological order 
of the phase. 

We emphasize that the expansion used here cannot be 
treated as a perturbative treatment in the adiabatic limit where 
the toric code coupling constants are turned on slowly because 
in that case X < h for short times. We know that in the case 
of adiabatic evolution, we approximately remain in equilib- 
rium, i.e. in the ground state manifold, and topological order 
emerges after a time of order system size^ 



V. CONCLUSIONS 

We studied the quench dynamics of the entanglement en- 
tropy in the Kitaev toric code model. In the case of the inte- 
grable quench we found the entropy exactly as a function of 
time. The leading order contribution was found to scale as 
the subsystem perimeter with an oscillatory behavior that gets 
suppressed in the presence of some disorder in the coupling 
constants leading to the saturation of the entanglement en- 
tropy to disorder-independent values respecting the area law. 

We obtained results for the generated entanglement entropy 
as a function of an initial paramagnetic state with arbitrary 
magnetization. We studied the effect of the sharp corners in 
the subsystem boundary and found that in analogy with the 
static entanglement discussed in the literature, they make sub- 
leading contributions to the dynamically generated entangle- 
ment entropy as well. We showed that the topological entan- 
glement entropy remains zero in this quantum quench. The 
results of the integrable quench as well as the analysis of the 
integrability-breaking perturbations strongly suggest that the 
non-equilibrium state after a quench into the topological phase 
does not develop any of the topological structure present in the 
equilibrium phase. 

We used the concept of the entanglement entropy in the in- 
teraction picture which simplifies perturbative approaches to 
the calculation of the entanglement entropy. This allows us 
to tease apart in this particular problem the effects of the time 
evolution with an unperturbed Hamiltonian that cannot gen- 
erate either entanglement with volume scaling or topological 
entanglement and focus only the part of the evolution that can 
possibly lead to volume scaling or topological entanglement. 
By studying the effects of an integrability-breaking perturba- 
tion, we found a mechanism where volume scaling of the gen- 
erated entanglement entropy can emerge at high orders in per- 
turbation theory as expected of thermalized states. The notion 
of thermalization after a quench naturally ties in with topolog- 
ical order through the fragility of topological order to thermal 
fluctuations. 
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Appendix A: Derivation of entropy for y direction initial 
magnetization 

We explicitly calculate in this appendix the contributions 
to the entanglement entropy from regular boundary stars and 
plaquettes as well as the convex and concave corners when the 
initial magnetization is in the y direction. 

First let us find the number of regular boundary stars iSr) 
and plaquettes (P^) in terms of the subsystem perimeter and 
the number of convex and concave corners. We define the 
perimeter P as the total number of boundary plaquettes or 
boundary stars (they are equal). We then have 



Pr=P- C-, 



2C< 



Sr=P- 2C^ 



c^ 



Once again we use Eq. (|8]l which generally holds for any 
initial state jfi). The only intermediate states (in the a^ ba- 
sis) which contribute to the tr(p^) are the ones connected to 
|J7,y) by some star or plaquette flips. We label these states by 
variable = ±1 and -0 = ±1 living on the boundary plaque- 
ttes and stars respectively. A negative Q (ip) indicates that the 
corresponding plaquettes (star) is flipped with respect to \^)y 



\al3) = \{ep^,^s 



9p,>^^,V.<,0iJ) 



(Al) 



where Pr, Sr, p> and s< label regular plaquettes, regular 
stars, convex corner plaquettes and concave corner stars re- 
spectively. For j — 1,2, {M } ({9i }) are variables liv- 
ing on the two stars (plaquettes) corresponding to the convex 
(concave) corners. Notice unlike regular star variables the two 
stars at a convex corner are labeled by the corner plaquette and 
similarly the two plaquettes at a concave corner are labeled by 
a corner star since a convex (concave) corner uniquely corre- 
sponds to a corner plaquette (star). 

Now suppose two such states have the same (3 spins. All 
regular and concave corner plaquette and star flips must be 
the same for these two states. As seen in Fig. |2] however, at a 
convex corner the 9 or i/j variables for these two states could 
be different if both the convex corner plaquette and the corre- 
sponding two corner stars simultaneously have opposite flips. 
Then for two states |a/3) and \a' (3) (which can be labeled as 
Eq. ( lAlb because they contribute tr{p^)) all their 6 and ip 
variables must be the same except for the ones correspond- 
ing to convex corners, namely {6'p^, V'^ } and {dp^,ipl }' 
which follow the relation 

{^p> ' ^p> y = {ip> &p> ' 7p> ^p^ } 

with auxiliary variable ^p^ = ±1. Similarly we can introduce 
variables ps^ — ±1 for concave corner and for two states 
I a/3) and |a/3') we must have 



Introducing a replica index I on 7' and p', we can write 

|ai/3i) = |{0p^,7A,^,^p,,V^^,V.„,0i}) 

|a2/3i) = |{0p,,^,,,7p^0p>,7p>V'^>,V'.<,^iJ> 



Let us for the moment focus on the Renyi entropy — tr(p^) 
for simplicity. The calculation for tr(/9^) will be very similar. 
First we can see from the form of the the states above that the 
sign factors from Bp acting on states in the a^ basis have no 
effect. A Bp plaquette flip can only give a different sign for 
the matrix elements involving |ai/3i) and |q!2/3i) at a corner. 
But then the plaquette flip at that corner would also give a 
different sign for the matrix elements involving \a1li2) and 
\a2P2) and the overall sign for the product of the four matrix 
elements appearing in tr(/9^) will be positive. 

The calculation is then similar to the simpler cases with the 
initial magnetization in the x and z directions except for each 
convex or concave corner there is an additional auxiliary vari- 
able to sum over Again, each matrix element can be written 
as a product similar to Eq. ( fTTT ). By simply bringing together 
terms we can write tr(p^) in the following form: 



tr(py 



T.\{U.\{fs.\[fp>\{fs^ (A2) 



p> 



where the summation is over all 6, ip, 7 and p variables for 
both replicas. We can take the sum Eq. lA2l inside the products 
by switching their order and perform the sum for individual 
terms. This gives the factorized form below 

tr(pi) = n ^P'- n p^'- n pp> n ^«< ■ ^^^^ 



where we will compute the Fs below. Assuming no disorder 
and taking the log of both sides, we immediately obtain our 
main result for the entropy 

S= (P-C>-2C<)lnFp^ 

+ {P- 2C> - C< ) In F,„ + C> In Fp^ + C< In F^^ 

which conforms to the general expression Eq. ( fTSl l. 

For regular stars and plaquttes, the terms appearing in the 
above Eq. ( IA2I ) are identical to what we obtained when the 
initial spins were in the x or z direction. For example 



fp.= 



f^ + ^P,- 2^ . 

I — - — COS Xp^t + 



1- 



— sin XpJ 



Let us now define the following function for a more com- 
pact notation 



li^,x) 



l + § , , 1-^^. 2 



■ cos X 



sm X. 



We can then write 

/p. = [g(V, Ap^i)]' , /,, = [gii:s^,p,J)f . (A4) 
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At the corners, we have a more comphcated structure and 
for example, fp^ , in addition to 9p^ , also depends on 7^ and 
ipijj for j — 1,2. Notice that for the Renyi entropy a^ = ai 
and therefore 7^ =1. We can explicitly write 

n 5'(^p> ' ^^p> ^)9(^P> ^i> ' ^^i> *) ■ 

It can be easily shown that 

g{^,x) = |sina;cosx|e'''"l™*^L 
Using the above identity we can write 

fp> = (gSin2Apj.t sin2^p_^i sin2/ip^t)^ x 

g(l+7p>)(ep>ln|cotAp>t|+Ej'AJ>ln|cotpJ^t|)^ 

A similar expression with 9 ^ ip, X ^ fi and ^ ^ p gives 
fs^ ■ We can now perform the summations and we obtain 



^p. = H /p. = (cosAp^t)^" + (sin Ap,_i)2 



(A5) 



with n = 2 for the Renyi entropy and similarly for F^^ with 
A -(-^ /i. To calculate 

^p> ~ 2^ ■'P> ' 

we first perform the sum over ^s and ^s and then over 7s and 
we obtain after some algebra 

Fp^ = ^ (— sin^ Xp^t sin^ /i^^i sin^ /ip>* ^^6) 
fc=i ^'^ 



„2 ,,2 .\" 



+ ^feCos Ap^t COS ^p^t COS iip^t 

with n = 2 for the Renyi entropy and ^1 = 1, ^2 = tan^ Ap-^ t, 
^3 = tan^ /ip t and ^4 = tan^ /i^ i. We can similarly find 
Fs^ for concave corners by switching X -^ ji. The general- 
ization to n > 2 is similar and leads to the appearance of a 
different n in Eq. (IA6I 1 and Eq. 



Appendix B: Transfer matrix method for an arbitrary Initial 
magnetization 

In this appendix, we discuss the transfer matrix method 
used in calculating the Renyi entropy for an arbitrary initial 
magnetization in the positive n direction. The initial state can 
be written as 

l^> = (c+|+>+c-|-»^"' 
where |+) and |— ) are the eigenstates of cr^. 



Let us begin by writing an explicit equation for tT{p']^. Us- 
ing n = 2 in Eq. ([8]l we can write 



Hpa 



^(ai/3i|C/+(7+|r!)(r!|C/7[/-|a2/3i) (Bl) 

{a2h\U+Ut\n){msU;\^iP2). 

Consider the first of the four matrix elements appeasing in the 
equation above, namely {aij3i\UpU'^\Vt). We represent the 
|ai/3i) spin configuration in a basis where boundary spins, i.e. 
spins belonging to a boundary star or plaquettes, are in the a^ 
basis and all other spins in the cr^ basis. In a schematic no- 
tation explained below, we can then write the matrix element 
as 



(ai/3i|t/+C/+|f])=^[, 



i^t ^ a'/3aa 



{6} 



(B2) 



(cosAi)'^M*sinAi)^^ c+'c 



where in cos Xt and sin At we can more generally make the 
substitution At — > /^ X(t')dt' . 

The notation in the above equation needs some explana- 
tion. As before we have introduced Ising variables 6p on the 
P boundary plaquettes. For a given configuration of these 
variables {0}, p^ and p^ are the number of flipped and un- 
flipped boundary plaquettes respectively and b^ and b^ are the 
number of up and down boundary spins in the state obtained 
from |ai/3i) after flipping the plaquettes represented by {0}. 
We have introduced a new notation where a denotes only the 
boundary spins in subsystem A that belong to a boundary pla- 
quette and the other boundary spins in A, which only belong to 
a boundary star, are denoted by a' . Similarly boundary spins 
in subsystem B that belong to a boundary star are labeled /3 
and the rest /3' as seen in Fig.|5] Assuming that the total num- 
ber of boundary spins is B, we can relate p^,^ and 6^ 4^ to the 
spins variables above as follows 



P^+Pi^ P, Pt - -P4. = X! ' 



and 



&t + ^i 



B 



5:.'+5:^a+5:^/3'+5: 



In the sums above each a and /?' spin is multiplied by 
the Ising variable 9 of the boundary plaquette it belongs 
to. Each (3 spin belongs to two boundary plaquettes and is 
therefore multiplied by the two corresponding 9 variables. 
In the notation explained above, we now write the term 
(ai /3i I [/+ [/+ 1 17) as follows 



(a 1^1 1 [/„+[/+ 1 17) = {c+c-)- (i sin At cos At)- e 



T pifJ-tJloi' flaa 



X A^ ^ "' V v4^(5: ea+E 0/3'+E ope) ^i E i 

{6} 



where 



H 



-i cot At, 



A: 
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FIG. 5. The notation used in Eq. lB2l The boundary spins in subsys- 
tem A are divided into a' spins that do not belong to any boundary 
plaquette and a spins that do belong to a boundary plaquette. Sim- 
ilarly, the boundary spins in subsystem B are decided into /3' spins 
that do not belong to any boundary star and /3 spins that do belong to 
a boundary star. We have represented the a' spins by a black circle 
on a dashed black bond, the a spins by a black circle on a solid black 
bond, the /?' spins by a blue square on a dashed blue bond and the /? 
spins by a blue square on a solid blue bond. 



Like the one-dimensional classical Ising model in a mag- 
netic field, we can now integrate out the 9 Ising variables using 
the transfer matrix method. The terms containing 9 are pro- 
moted to 2 X 2 transfer matrices where each transfer matrix 
^2x2 is between two adjacent boundary plaquettes and 

{ail3i\U+U+\n) = (c+c-)^ {i sin Xt cos Xt)^tT(j[T2^2 

(B3) 
The other degrees of freedom can be systematically integrated 
out by repeatedly using the following trace identity 



tT{AiA2 ■ ■ ■ )tT{BiB2 



\ ^ tT[{Ai (g> Bi) {A2 (g) B2) ■ ■ ■] . 

(B4) 
Repeated use of the identity Eq. (IB4b leads to enlarged trans- 
fer matrices as discussed below. Tracing out all degrees of 
freedom will lead to 64 x 64 transfer matrices the trace of 
whose product gives the Renyi entropy. 

Notice that at concave and convex corners, we get a differ- 
ent 2x2 transfer matrix. The 2x2 transfer matrices at the 
corners will be given at the end of this appendix for complete- 
ness. We discuss below the procedure for obtaining the final 
64 X 64 transfer matrices from the 2x2 ones along a smooth 
boundary. At the sharp corners, we can obtain the 64 x 64 
transfer matrices from the corresponding convex and concave 
corner 2x2 matrices essentially in the same way. 

Let us now explicitly write T2X2 in Eq. (IB3l l along a smooth 
boundary. For a bond /3* between two adjacent boundary pla- 
quettes as show in the top left corner of Fig.|6] we have 



To 



A2" 



exp 



t^ta p a a 



i+i 



H-iA-i<^0"+'^'^ 



AhP' A-hf^' 
A-W A^' 




nL>. 








Ql+l 



^R>. 



FIG. 6. Top left: the transfer matrix 72x2 along a smooth boundary 
and the corresponding spin variables. Top right: the transfer matrix 
r^x'2 for entering a convex comer. Bottom left: the transfer matrix 
2^2 x'2 for exiting a convex corner Bottom right: the transfer matrix 
2^2'x 2 foi' 2 concave corner. 



The next step is to consider 

i?i^2 = Y.{aiPi\u+ut\n){n\u-u-\a2^i) 
Pi 



We can now use the trace identity Eq. IB4I to bring together 
terms with the same /3 and /3' and trace out these degrees of 
freedom to get 



Ri,2 = (|c+c"|)^(sinAtcosAt)^tr (n^4x4) 
where the 4x4 transfer matrix above is given by 



Tix4 = exp i 13 fit (a'lalal'^^ - ci\oL2Ci'2^\ 
with the following ti matrices 

^ i/5 \ ( H*^ 



t\ t2 ts t4 



U 



t2 = 



H-2 

A^"! 







A- 



^ f A-2 0' 

I3'=±l ^ 




A* 2 



0' 



/3=±1 



AhP A-'2P 






Using the trace identity above to bring together the degrees 
of freedom via the tensor product of matrices is the main strat- 
egy for calculating the Renyi entropy. These calculations look 
tedious but are greatly simplified using Mathematica. The 
next step is to bring together the same degrees of freedom 
in i?i, 2^2,1 using the trace identity which results in 16 x 16 
matrices that only depend on a and a' degrees of freedom for 
both replicas. Notice that so far we have integrated out all (3 
and /?' degrees of freedom (for both replicas). 
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Now since the a' degrees of freedom do not interact we 
can do the sum directly. The a degrees of freedom have an 
Ising type interaction so the transfer matrix method can be 
applied again promoting functions of a (which appear in 16 x 
16 matrices) to 4 x 4 matrices (because there are 2 replicas). 
This operation can be done very simply in Mathematica. At 
the end of the day the Renyi entropy will be given by the trace 
of the product of P (the number of boundary plaquettes) 64 x 
64 transfer matrices. 

If the boundary of the subsystem is largely smooth, i.e. only 
has a few sharp corners, the different transfer matrices at the 
corners will give an order-one contribution to the Renyi en- 
tropy and the leading order area law term is 

^2(pA)«-PlogA 

where A is the largest eigenvalue of the smooth boundary 64 x 
64 transfer matrix constructed above. 

An exact calculation of the entanglement entropy is needed 
for evaluating the topological entropy. Therefore corrections 
from concave and convex corners must be taken into ac- 
count by using different transfer matrices at the comers in the 
boundary. There is a small complication due to the the fact 
that a convex (concave) corner plaquette has an additional a 
(/3) spin. Let us review our method: we first consider two ma- 

I 



trix elements each given by the trace of a product of 2 x 2 
matrices, we then bring together the same degrees of freedom 
and explicitly sum over /3s to get the 4x4 transfer matrices. 
Now having two /3s to sum over does not change anything ex- 
cept we have to do a few more summations to get r4x4. With 
our method of first integrating out the B spins and then the A 
spins, the additional a spin at a convex corner gives a jump in 
the indices of the transfer matrices at the end of the day but 
this can be easily accounted for by inserting a matrix Cy — 1 
between the transfer matrices and calculating the trace of their 
product as usual. 

Let us now write the 2x2 matrices for the corners. First 
consider the transfer matrix T^^^ for entering a convex corner 
as seen in the top right corner of Fig.|6] Each transfer matrix 
is written for a /3 bond and the a spins before (ctb) and after 
Qfa that bond. As mentioned above the discontinuity in the a 
index (aa for a transfer matrix is not equal to ab for the next 
transfer matrix) can be taken care of by inserting the matrix 
li6xi6'8'C4x4(4x4 because we have two replicas) at the end 
of the day between the 64 x 64 transfer matrices. For exiting 
a convex corner plaquette, with the notation of the bottom left 
corner of Fig.|6] we can write a 2 x 2 transfer matrix T^^-^ and 
similarly a transfer matrix T'j'^j for ^ concave corner (bottom 
right corner of Fig. |6]l. In terms of the spin variables shown 
on the corresponding corner of Fig.|6] these matrices are given 
by 



r2^>2 = A^""exp[zAiia"/3^aj,al 
T^^2 = ^^"" exp \i^it a'"/3'alai 



i/^A^('3"+"b) \ / AiP' A--2P' 

H-^2A-W'+<) ] \ A-hP' A^f' 



H^A^"^ 

H~^A-h< 



A-W AhP' 



ij5A3(ft+'3"+"b) 







iJ-iA-5(ft+/3"+«J,) 




A5"a 



A 1 fl' + l t 1 fl' + l 

A2P A^2P 

A2P A^2P 



A-^" 



r 



Notice that for a concave corner, the transfer matrix is writ- 
ten between two neighboring boundary plaquettes which in 
this case do not have a common edge because the plaquette in 
the middle lies entirely in subsystem B and is therefore not a 
boundary plaquette. The procedure for going from the initial 
2x2 transfer matrix to the final 64 x 64 is analogous to the 
smooth boundary case discussed above and a detailed discus- 
sion would not be illuminating. We have explicitly verified 
using the transfer matrix method described in this appendix 
that the order-one correction to the area law is entirely from 
the convex and concave corners and the topological entangle- 
ment entropy remains identically zero. 



with the spins inside (but not on the boundary of) the blue 
regions belonging to subsystem B and all other spins to A. 
For partition 1 for example, the outer and inner boundaries of 
the subsystem B are squares of side li and £2 respectively. 
The whole system is on an i? x £ torus. In terms of these di- 
mensions, we can explicitly find the quantities rip^, rips, ndA 
and n^B for the four different partitions and we observe that 
the value of each quantity for partitions 11 and 111, which are 
obviously equal, is the average of those for partitions 1 and 
IV. This immediately implies that the topological entangle- 
ment entropy identically vanishes at second order in h. For 
completeness we include the results of this simple geometry 



Appendix C: Renyi entropy at order /i"' 



We stated in section|IV]that Eq. ( |22] | gives zero topological 
entanglement at order h^. Consider the partitions of Fig. [3] 
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FIG. 7. The system for o(/i^) calculation. The dashed green bonds 
belong to subsystem B and the solid black bonds to the subsystem 
A. 



exercise for partitions I and IV below 



lA=^^-il+e\ 






{£i - 2'f -ll- 4£2 



n\B - 2 (4 - 5^1 + ^? - 5^2 - ^ 

^■pA 
"pB 



(^i-2)(£i-^2-4) 

2^2 _ 2^2 ^£2 +4 (2^2 
n^^ = 24 + 2£2_|_5^2_£^(^5. 



^3) 

2^2) 



where we have assumed £2 ~ ^1 > 4. 

Next, we discuss the fourth order term. The combinatorics 
get more complicated at higher orders but we can calculate 
the terms in the expansion by direct enumeration. This gives 
semi-analytical results with the time-dependence entering an- 
alytically through the multiple integrals depending on func- 
tions c, p and d and the geometry dependence computed for 
a system of given dimensions by computer enumeration. We 
have performed these calculation for several systems. Here 
we present the details for an 8 x 8 system shown in Fig.|7] 

If tr(p^) = 1 + /ft.2 + y/i4 + . . . ^ we have 

Let us assume we have coefficients fi and gi for partition 
i — 1 • • • 3. We can then write the fourth order term in the 
topological Renyi entropy as 



^topo 



o(M) 



h^ [252 



51-53 



Wf2 



fi~fl)] 



To find the fourth order term in S'topo, in addition to the 
/ coefficients that we found before, we need to calculate 
the coefficients gi, 32 and 173 for the different partitions. 
Looking at the structure of the Eq. ( IB lb . we observe that 
at fourth order we need to insert four Vs into the four bins 
(■■■)(''■)(■'')(■'■) which give a total of 35 terms of five 
different categories. There is one term of type (1, 1, 1, 1) 
(with one V in each (•••)) and 12 (2, 1, 0, 0), 6 (2, 2, 0, 0), 



12 (3,1,0,0) and 4 (4,0,0,0) terms. First consider the 
(1,1,1,1) term. We define a matrix A such that 



•^a,/: 



and obtain 



xp\V\Q) 



{AA^y 



5(1,1,14) = t'' 

Now since one V term acting on |0) can only change the state 
by flipping either one plaquette or one domino, the \a(i) states 
with a non- vanishing Aa.p can be enumerated and the matrix 
A can be constructed. The elements of the matrix A contain 
integrals (c), {p) and [d) depending on which term in V con- 
tributes to the matrix element. We then obtain for example the 
following expression for partition 1 of the 8x8 system 

51,(1,1,1,1) = 8 [33554432(c)(c)(c)(c) + 299(d)(d)(d)(d) 
+ 6336(d) (d)(p)(p) + 34176(p)(p)(p)(p) 
+ 196608(c)(c)(<i)(d) + 2359296(c)(c)(p)(p)]. 

We can similarly calculate 5(1,1.1,1) for other partitions. The 
next two categories are (2, 1, 1, 0) and (2, 2, 0, 0). The expres- 
sion for 5(2,1,1,0) and 5(2,2,0,0) are given below. 

5(2,1,1,0) = 2 ( - mVV\0){Q\V^\aP){aO\V\0) 
- (a0|yy|0)(0|y1'|a/3)(0/3|y|0) 

+ (a/3|yy|o)(o|yi'|o/3)(o|i^Vo) + c.c.) 



and 



5(2,2,0,0) 



{{<dP\VV\Q){Q\V^V^\Qp) 4- 

(o|yy|o)(o|yy|o) + c.c. ). 



:tO|yi/|0)(0|Wl/VO) 



Computing these expressions in terms of the functions c, d 
and / with direct enumeration yields the following expression 
for the topological Renyi entropy. 

^topoL(;,.) - -128 {I6(d)(p) [{p,d) + (rf,p)] 

+ {df-A{d,df-lQ[{p,d) + {d,p)f 
-48(p)4 + 256(p)2(p,p) - 320(p,p)2}. 

Notice that the last category 5(4,0,0,0) i^ independent of 
the partitions and does not contribute to the topological en- 
tanglement entropy at fourth order in h. We can show that 
the expression in Eq. dCll l above identically vanishes using 
[p.p) — (p)'^/2, {d,d) — {d)^/2 and the following integral 
identity 

{p,d) + id,p) = id){p). 

The fourth order calculation also confirms the notion that 
a growing ring-like (with the expansion order in h) region 
around the system boundary contributes to the entanglement 
entropy. 
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